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Abstract 



^^ ■ Finding the eigenvalues of a Sturm-Liouville problem can be a computationally 

(~| ! challenging task, especially when a large set of eigenvalues is computed, or just 

"j^ \ when particularly large eigenvalues are sought. This is a consequence of the highly 

oscillatory behaviour of the solutions corresponding to high eigenvalues, which forces 
a naive integrator to take increasingly smaller steps. We will discuss some techniques 
that yield uniform approximation over the whole eigenvalue spectrum and can take 
^ ' large steps even for high eigenvalues. In particular, we will focus on methods based 

|0 . on coefficient approximation which replace the coefficient functions of the Sturm- 

vQ \ Liouville problem by simpler approximations and then solve the approximating 

CNl ' problem. The use of (modified) Magnus or Neumann integrators allows to extend 

t;;j- . the coefficient approximation idea to higher order methods. 

o 

oo 
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1 Introduction 



The classical Sturm-Liouville problem (SLP) consists of a linear second-order 
ordinary differential equation written in formally self-adjoint form 



- {p{x)y')' + q{x)y = Xw{x)y, 
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defined over an interval a < x < b with appropriate boundary conditions at 
a and b. An eigenvalue is a value of A for which (1) has a nontrivial solution 
y subject to the boundary conditions, and the solution, unique up to scalar 
multiples, is the associated eigenfunction. 

The regular theory which dates back to Sturm (1809-1882) and Liouville (1803- 
1855) assumes that the coefficient functions are well-behaved, say p{x),q{x) 
and w{x) are piecewise continuous with p and w strictly positive, on a bounded 
closed interval [a, b], and that regular boundary conditions are imposed, namely 

aiy{a) + a2p{a)y'{a) = 0, biy{a) + b2p{a)y'{a) = 0, (2) 

where ai, a2 are not both zero, nor are bi, 62- Then there is an infinite sequence 
of eigenvalues 

Ao < Ai < A2 < . . . 

and eigenfunctions 

yo{x),yi{x),y2{x),..., 

such that yk{x) has just k zeros on the open interval {a,b), and such that 
distinct eigenfunctions are orhogonal with respect to the weight w{x): 

b 
yi{x)yj{x)w{x)dx = 0, i^ j. (3) 



The mathematical theory of SLPs is immense (see e.g. [1]) but they are not 
just objects of interest to mathematicians alone. Since the early 19th cen- 
tury SLPs have been ubiquitous in applied mathematics since they are the 
one-dimensional models of a large number of important physical processes in 
fields such as acoustics, geophysics, waveguide theory, hydrodynamic stability, 
neutron transport, .... They arise in the analysis of such processes in more 
than one dimension by the method of Separation of Variables. Another reason 
why SLPs are of vital interest to physicists is that Schrodinger's equation in 
one dimension is of Sturm-Liouville form. 

Some eigenvalue problems have explicit solutions, and are therefore impor- 
tant in the analytical investigation of different physical models. However most 
eigenvalue problems are not solvable, and computationally efficient approxima- 
tion techniques are of great applicability. The numerical solution of (regular) 
Sturm-Liouville problems is not trivial. The choice of numerical method for 
efficiently approximating a sequence of eigenvalues of the SLP depends on the 
desired accuracy of the estimates and also upon the number of eigenvalues 
required. General ODE boundary-value software can solve SLPs, but ineffi- 
ciently. The challenges are to do this more cheaply, especially when long runs 
of higher-order eigenvalues are required. 
In fact, many classical methods involve the approximation of the corresponding 



eigenfunctions by piecewise polynomials and are thus inefficient for the com- 
putation of higher eigenvalues which have severely oscillatory eigenfunctions. 
There have e.g. been many developments in the basic approach of reduction 
to a matrix eigenproblem using finite differences and finite elements. Excellent 
surveys of such matrix methods are given in [2,3,4]. Matrix methods can only 
be used to approximate the ffist few eigenvalues. The error for even moderately 
large k is considerable unless the dimension of the associated matrix is very 
large. In this paper we will concentrate on a different class of methods which 
are based on shooting-type algorithms. These methods perform much better 
than the matrix methods for singular (or nearly singular) problems, for the 
computation of eigenfunctions, and even for the highly accurate computation 
of the ffist eigenvalues. We will discuss in particular some important contri- 
butions to the efficient and accurate computation of the higher eigenvalues of 
SLPs. 



2 Shooting methods 



Shooting methods are based on the reduction of the boundary value problem 
(l)-(2) to the solution of an initial value problem. The differential equation is 
solved as an initial value problem over the range [a, b] for a succession of trial 
values of A which are adjusted till the boundary conditions at both ends can be 
satisfied at once, at which point we have an eigenvalue. The simplest technique 
is to 'shoot' from a to b. This means that one chooses initial conditions which 
satisfy the boundary condition (2) in a: 

yda) = -02, p(a)|/^(a) = ai 

The boundary condition at b determines 'target' values; if the value of y 
matches the target, we have found an eigenvalue. In fact, the eigenvalues are 
determined, using some iterative technique, as the solution of 

0(A) = %l(6, A) + M%l(^ A) = 0. 

Another option is to shoot from the two ends to some matching point a < 
Xm < b. In this case also a right-hand solution yji is defined satisfying the 
conditions 

ydb) = -62, p{b)y'L{b) = bi. 

The two solutions yi and y^ are arbitrarily normalised, so their values can 
always be made to agree at the matching point by renormalising them. How- 
ever, the criterion for a trial value of A to be an eigenvalue is that the ffist 
derivatives should match, as well as te values. The mismatch function (also 
called miss-distance) is thus given by 

0(A) = yiixm, \)p{xm)yR{xm, A) - ynix^, \)p{xm)y'L{x^, A). 



Numerically some choices of Xm may make it more difficult to compute 0(A) 
than do others (see [5]). Generally it is a good idea to take the matching point 
in the interior of the interval, away from singular endpoints. 

Thus the procedure for finding the numerical value of an eigenvalue, consists in 
evaluating the mismatch function 0(A), numerically, and then through a finite 
series of iterations finding the value of A such that 0(A) = to the required 
accuracy. The usual iterative methods for finding the roots of a function may 
be employed here to find the zeros of 0(A). One problem associated with 
this approach is that the function 0(A) does not give any way of determining 
the index of the eigenvalue once it has been found. Thus we have no way of 
knowing which eigenvalue we have found when 0(A) = 0. Likewise, in order 
to converge on a specific eigenf unction, one has to enhance the algorithm, for 
instance by counting the zeros of the solution as part of the integration for 
each trial A value. 

To solve this problem, the polar coordinate substitution, known as the Priifer 
transformation, is used. This transformation makes it possible to specify the 
eigenvalue to be computed. The (scaled) Priifer transformation is defined by 
the equations 

y = S-^'^p sin e, py' = S^'^p cos 9, (4) 

where 5 is a strictly positive 'scaling function' chosen to give good numerical 
behaviour. In [5] it shown that the resulting differential equations for p and 9 
are then of the form 



9 = — cos'' 9 + ^ ^^-^ sm^ 6* + — sm 6* cos 9, (5) 




, sin 2^ -—COS 2^. (6) 

p S J S 

The regular boundary conditions (2) at a and b define the conditions for 9 

9{a) = a, 9{b) = (3, 

where 

S{a)a2 S{b)b2 

tana = , tanp — 



Ol bi 

These equations only determine a and /9 up to a multiple of vr. As stated by the 
following theorem (proved in [5]), each (appropriate) choice of this multiple 
specifies in fact precisely one eigenvalue. 

Theorem 2.1 Consider the scaled Priifer equations of a regular Sturm- Liouville 
problem. Let the boundary values a and (3 satisfy the following normalization: 

ae[0,7r), /5e(0,7r]. (7) 



Then the kth eigenvalue is the value of A giving a solution of (5) satisfying 

e{a;X)=a, e{b; X) = (3 + kn. 

This leads to several useful numerical methods based on some form of the 
Priifer transformation. Priifer based shooting methods can be constructed 
where the counting of the zeros of y{x) needed to compute the specific eigen- 
value with a given index k is built in. One can for instance define a shooting 
method for the 6 equation: For any A, let Oi{x; A) and 6r{x; A) then be the 
solutions of (5) satisfying 

e^ia; \) = ae [0, tt), ^^(6; A) = /? G (0, n]. (8) 

The scaled Priifer mismatch function is then defined by 

(pi\)=9Lix^;X)-9Rix„,;\). (9) 

and the eigenvalue A^ is the unique value such that (p{\k) = kn. 

The SLEIGN code (and its successor SLEIGN2) from Sandia Laboratories 
[6,7] uses an (explicit) Runge-Kutta method to integrate the 6 equation. For 
certain problems, where a good scaling function S is heuristically found, the 
oscillations are removed and large steps can be taken. However there is no gen- 
eral method for finding a good scaling function and such shooting methods 
based on standard initial value libraries often suffer from stepsize restriction 
when solving for large eigenvalues. They also have some difficulties caused by 
stiffness of the 9 equation (5) in a 'barrier' region where {Xw — q)/p is large and 
negative. Instead of using a standard initial value library code, it is a better 
idea to combine a Priifer formulation with coefficient approximation, in which 
the coefficient functions are piecewisely approximated by low degree polyno- 
mials. Then the integrations may be performed analytically and stiffness is no 
longer a problem. 



3 Coefficient approximation 



An important class of methods for the numerical solution of Sturm-Liouville 
problems is based on coefficient approximation. The basic idea here is to re- 
place the coefficient functions p{x), q{x), w{x) of the Sturm-Liouville equation 
piecewisely by low degree polynomials so that the resulting equation can be 
solved analytically. The idea dates back at least to Gordon [8] and Canosa 
and De Oliveira [9] and was studied also by Ixaru [10], Paine and de Hoog 
[11] and Smooke [12]. But the standard reference is due to Pruess [13,14]. He 
examined the piecewise constant case and his strategy has been implemented 



by Pruess and Fulton in the code SLEDGE [15]. This so-called Pruess method 
replaces the SLP (1) by the approximating problem 

— {jyy'{x))' + qy{x) = \wy{x), x G (a, b) (10) 

where p, q, and id are piecewise constant (midpoint) approximations of the 
functions p, q, and w. The y{x) of the approximating problem (10) can then be 
integrated explicitly in terms of trigonometric and hyperbolic functions: Let p, 
q and w have constant values pi, qi, Wi in the ith interval (x^-i, x^), i = 1, . . . , n 
with step size hi = Xi — Xj_i: 

-iPiy'ix))' + qiy{x) = \wiy{x) 

the solution over [xj_i,Xi] is then advanced by the relation 




^{Zi) hir]o{Zi)\ ( y{xi_i) 
Zir]Q{Zi)/hi i{Zi) J ypiy'{xi-i)^ 

with Zi = hj{qi — \Wi)/pi and 

' sm(\z\^/^)/\z\^/^ a z <0. 



:iii 



1 if Z = , 



cos(|Z|i/2) ifz<0, 

aZ) = { Vo{Z) - 

cosh{Z^/^) if Z > , 

smh{Z^/^)/Z^/^ ifZ>0, 

(12) 
One can also propagate the solution from Xi to Xj_i, by taking the inverse of 
the transfer matrix in (11), which is just the result of replacing hi by — /ij in 
this matrix. 

This gives us a method for explicitly integrating {y,py') over the x range, 
and to use a shooting method. This is done e.g. in SLEDGE and combined 
with the ideas based on the Priifer substitution to be able to home in on a 
particular eigenvalue (see [15]). 

Pruess proved that if p, q and w are piecewise constant on a mesh of typical 
meshsize h and equal to p, q, and w at the mesh midpoints and if A^ is the 
kth eigenvalue of the approximating problem, then 

\^k — ^k\ ^ Ch k\Xk\ 

for all k and small enough h. Thus one would expect a higher eigenvalue to 
need more meshpoints to compute to a given relative tolerance than a lower 
eigenvalue. However, as mentioned in [5,16], there are two reasons why this 
is not seen in practice. Firstly, arguments show also that \Xk — ^k\ < Ch\Xk\ 
for large k and small enough h. Secondly, many problems occur in Liouville 



normal form (Schrodinger form) where p = w = 1 and for these there is an 
improved error bound 

\^k — ^k\ ^ Ch k \Xk\. 
Thus we can actually use larger h for large k for a given relative error. 

The Pruess-type methods have some important advantages. As already noted, 
Pruess methods are relatively unaffected by the stiffness/instability which can 
force a very small stepsize on a standard initial-value solver, and the accuracy 
is maintained (or even improved) as fc — > oo. Further advantages of the Pruess 
methods are that it allows a very simple interval truncation algorithm for 
singular problems (see section 5) and that unlike a method based on a standard 
initial- value solver, it is practical with the Pruess method to fix the mesh and 
evaluate the coefficient midpoint values once for all before the start of the 
shooting process. Since the overall shooting process consists of a number of 
integrations with different values of A the latter can give a big speed advantage. 

A drawback of the Pruess methods is the difficulty in obtaining higher order 
methods. It is usual to implement them using Richardson extrapolation. It 
is clear that the step sizes must be sufficiently small such that the error in- 
troduced by the approximation by piecewise constants is not too large. This 
means that for problems with strongly varying coefficient functions the num- 
ber of intervals in a mesh can be quite large. Some approaches have been 
suggested towards the realization of higher order methods based on coefficient 
approximation. These approaches can be classified as "modified integral series 
methods", which will be discussed next. 



4 Modified integral series methods 



We will consider two integral series which allow the natural extension of the 
Pruess-ideas to higher order methods: a Neumann series and a Magnus series. 
In fact, these integral series offer an easy way to approximate the coefficient 
functions of the SLP by higher order (piecewise) polynomials, giving more 
accurate results than the approximation by a piecewise constant. 



4-1 The Neumann and Magnus expansion 



There is an emerging family of numerical methods based on integral series 
representation of ODE solutions. Consider the linear differential equation 

N 



y' = A{x)y, y(0) = yoGM^. (13) 



The simplest integral series is obtained by applying Picard iteration [17] to 
obtain the fundamental solution of the matrix linear ODE 



y[x) 



rx rx rxi 

1+ A{xi)dxi + / A{xi) / A{x2)dx2dxi 
Jo Jo Jo 

rx rxi /■X2 

+ / A{xi) I A{x2) / A{x^)dxzdx2dxi + . . . 
Jo Jo Jo 



(14) 



yo 



This series is known as the Feynman-Dyson path ordered exponential in quan- 
tum mechanics, in mathematics it is known as the Neumann series or Peano 
series. 

The Magnus and Cayley expansions are two other examples. They are obtained 
by transforming Eq. (13) to the suitable Lie algebra and applying the Picard 
iteration to the transformed ODE. Details on both approaches can be found 
in [18]. The Cayley expansion is based on the Cayley transform while the 
Magnus expansion is based on the exponential map. The approach of Magnus 
[19] aims at writing the solution of Eq. (13) as 

y{x) = exp(fi(a;))yo 

where Q{x) is a suitable matrix. The Magnus expansion says that 



Q{x) 



A{xi)dxi 

2 Jo 

1 

+ 4.0 

1 

^12 Jo 



A{x2)dx2, A{xi) 



dxi 



A{x3)dx3, A{X2) 



dx2, A{xi) 



x\ 



A{x2)dx2, 



Xl 



A{x3)dx3, A{xi) 



dxi 

dxi + 



(15) 



where [■, ■] denotes the matrix commutator defined by [X, Y] = XY — YX . 

Numerical schemes based on the Magnus expansion received a lot of attention 
due to their preservation of Lie group symmetries (see [18,20] and references 
therein). The Neumann series does not respect Lie group structure but avoids 
the use of the matrix exponential. The use of Neumann series integrators has 
been proved successfull for certain large, highly oscillatory systems in [21]. 

Since the SLP can be written in the matrix form (13), both Neumann and 
Magnus schemes can be considered for the numerical solution of the SLP. The 
Sturm-Liouville equation (1) in matrix form reads 



y'(a;) = A{x)y{x) 



l/v{x) 

q{x) — \w{x) 



y[x) 



(16) 



with y-^ = {y{x),p{x)y'{x)). Moan [22] was the first to consider a Magnus 
series integrator for the SLP in the Schrodinger form y"{x) = {q{x) — \)y{x) 



or in matrix form 

( 



// 




y\x}= \y{x}. (17) 

\q{x) 

He applied the Magnus integrator directly to this problem. However poor 
approximations were obtained for the higher eigenvalues, as a result of the 
finite radius of convergence of the Magnus series [23]. When the solution of 
a linear system y' = A{x)y oscillates rapidly, modified schemes should be 
used, as recommended in [24,25,26]. Describing these modified schemes we 
will focus on the basic Schrodinger equation y"{x) = {q{x) — X)y{x), but 
the schemes can be extended to the more general Sturm-Liouville problem 
— {p{x)y'{x)y + q{x)y{x) = \w{x)y{x). 



4-2 Modified Neumann and Magnus schemes for the Schrodinger equation 



We consider the Sturm-Liouville problem in Schrodinger form eq. (17) which 
is a problem of the form 

y{xy = A{x,\)y{x), y{a)=yo, (18) 

where y = [y{x),y'{x)]'^. Note that the coefficient matrix is in s/(2),i.e. the 
matrix has a zero trace. 

Suppose that we have already computed yj_i ^ y{xi-i) and that we wish to 
advance the numerical solution to Xi = Xi-i + hi. We first compute a constant 
approximation q of the potential function q{x) 

q = — q{x)dx. (19) 

hi Jxi^i 

Next we change the frame of reference by letting 

y[x) = e(''"^"-^^^u(a; - Xi^i), Xi^i < x < Xi (20) 

where 

A{X)=\ ° M. (21) 

We treat u as our new unknown which itself obeys the linear differential 
equation 

u'((5) = 3(5, A)u((5), 5 e [0, k^, u(0) = y,_i (22) 

where 

B{S, A) = e"^^ {A{x,.i + 6) -a) e^^. (23) 



9 



The matrix B can be computed explicitly. With ^{Z) and tiq{Z) defined as in 
eq. (12) we can write B as 



B{6,\) = A,{6) 



^^°^^^^^ 2(A - q) 



6r]o{Z; 



(24) 



where Ag{S) = q — q{xi^i + 5) and Z^ = Z{j) = [q — A)7^. 

We have thus replaced one linear system by another. The new system (22) has 
one crucial advantage over (18): the entries of the matrix B are themselves 
rapidly oscillating functions (for X > q). This is not very helpful when (22) is 
solved by a classical method as e.g. a Runge-Kutta or multistep method. When 
the modified equation is however solved by an integral series method, repeated 
evaluation of integrals of B is required. This integration is a "smoothing" 
operator: the amplitude is decreased once the integrand is integrated. As a 
result the higher the oscillation, the faster the convergence of the integral series 
method and the faster the decay in local error. We can refer to [24,25,26] for 
numerical results confirming the success of this approach for highly oscillatory 
ODEs. 

Over each interval [xi_i,a;j] an integral series is applied on the transformed 
equation u'{6) = B{6)u{6). This requires the truncation of the integral series 
and the replacement of integrals by quadrature (see next section). The solution 
y in X = Xj is then obtained from y(A,Xj) = e^^'^u^hi). Note that e^^^ is the 
known solution of the system with constant potential 

, /.A / az,) hMZk) 

expm = \ , Zh = Z{hi) (25) 

^K{q-\) 0] [Zj,r]o{Z,,)/h, ^Z^) 

and thus the same as the transfer matrix in (11). 

The first option we consider is the use of a Neumann scheme. Application 
of the Neumann series integrator to the modified equation u'{6) = B{S)u{6) 
gives 

u(hi) = 1+ B(x)dx + / / B(xi)B(x2)dx2dxi 

I Jo Jo Jo ^26) 

+ / B{xi)B{x2)B{x3)dx3dx2dxi + . . . yj_i 

Jo Jo Jo 

When only the first term in the Neumann series is retained, one has u(/ij) = 
yi_i and with y{xi) = e'^^^u^hi) this is exactly the second-order Pruess method 
given by eq. (11). Higher order methods are obtained by including more Neu- 
mann terms. In [261 it was shown that in fact each extra Neumann term can 



10 



be seen as a correction term in a Piecewise constant Perturbation Method 
(PPM) of Ixaru and co-workers. The PPMs use a perturbation technique (the 
first systematic description of this technique is due to Ixaru [27]) to construct 
some correction terms which are added to the known solution of the approxi- 
mating problem y" = {q — \)y with a piecewise constant potential q. The more 
correction terms included the higher the order of the algorithm. In [28,29] the 
PPM algorithm is described and applied on regular Schrodinger problems. 
The PPMs were also extended to general Sturm-Liouville problems using the 
Liouville transform and formed the basis of the Fortran code SLCPM12 [30] 
and the graphical Matlab software package Matslise [31]. 
To approximate the integrals in (26) quadrature must be used which can deal 
adequately with the oscillatory entries of the matrix function B. In section 
4.3 a Filon-type quadrature rule will be discussed which is very similar to the 
procedure used in the description of high order PPMs in [28,32]. There the 
potential function q is replaced by a piecewise polynomial, which makes the in- 
tegrals in (26) analytically solvable. The degree of the (piecewise) polynomial 
can be taken sufficiently large such that this approximation has no influence 
on the accuracy of the method. The order of the method then only depends 
on the number of terms retained in the Neumann expansion (i.e. PPM cor- 
rection terms). As mentioned, including only the first Neumann term gives 
us a method of order two. Including also the second Neumann term (the first 
integral) leads to a method of order four, adding the third term (the double 
integral) results in an eighth order method and adding the fourth term gives 
us a method of order ten (see [21]). 

Another option is to apply a Magnus method to the modified equation (22). 
The Magnus expansion is then 



a{5) = ai{6) + ^2(5) + a^iS) + ai{6) + ..., 



(27) 



where 



.5 

c"i(^) = / B{x)dx, 
Jo 

0-2(5) = -- / / [B{x2),B{xi)]dx2dxi, 
2 Jo Jo 



^3(5) 



1 

V2Jo 
1 r^ 

I Jo 



X\ 



B{x2)dX2, 



Xl 



B{x2)dx2, B{xi) 



X2 



B{x3)dx3, B{x2) 



dx2, B{xi) 



dxi, 
dxi, 



f28) 



and u(S) = e^'^^^Yi^i, 6 > 0. Thus, to compute y^ = e'^^e'^'^^^yj^i with h = hi, 
we need to approximate cr{h) by truncating the expansion (27) and replacing 
integrals by quadrature (see 4.3). As shown in [26], truncating all but the 
first integral leads to a fourth order method, while including also (72 gives us a 
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scheme of order eight. Having approximated cr{h), its 2 x 2 matrix exponential 
must be computed. We note that ct(/i) is always a two by two matrix with zero 
trace. For such matrices the following is true: 

^{uj) + ar]o{uj) br]o{uj) \ 

I cu = a^ + bc. (29) 

cr]o{uj) ^{lo) - ar]o{uj) J 
Here a, b, c, u are functions of x and E. 

4-3 Quadrature of the (multivariate) integrals 



Practical implementation of both the Neumann and Magnus series requires 
the replacement of multivariate integrals by quadrature. Although multivari- 
ate quadrature is usually considered a hard problem, it is possible to imple- 
ment Neumann and Magnus expansions with surprisingly cheap and effective 
quadrature. Moreover when a Filon-type quadrature method is used, even the 
highly oscillating integrals, which appear when X ^ q, are approximated to 
a suitable precision in a small number of function evaluations per step. Filon 
quadrature has been analysed extensively in [33]. 

For a Neumann integrator as well as for a Magnus integrator, the univariate 
(modified) integral /q ' B{6)d6 needs to be approximated. A Filon-type rule is 
used. Here this means that Ag(5) in (24) is replaced by a polynomial, i.e. by 
the Lagrange polynomial 

£^^i5) = J2\icihM5) (30) 

1=1 

where ii is the /th cardinal polynomial of Lagrangian interpolation and Ci, 
C2, ... ,Ci, are distinct quadrature nodes. The resulting integrals can then be 
solved analytically. For each entry in the univariate integral a scheme of the 
following form results 

V 

hiJ2ki^)^qicihi), u = q-\. 
1=1 

When no further Neumann or Magnus terms are retained in the algorithm, the 
truncated Neumann or Magnus scheme is of order four and it is then sufficient 
to have v = 2 Legendre quadrature nodes (or z/ = 3 Lobatto nodes). This 
means thus that in this case A^ is approximated by a linear polynomial. 
For schemes of order eight, the double integral must be included and z/ = 4 
Legendre nodes should be used. As for the univariate integral, the double 
integral is computed by replacing A^ by the polynomial £Aq and solving the 
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resulting integrals analytically (using a symbolic software package). Each entry 
in the double integral is then approximated by an expression of the form 

V V 

^i 51 XI bk,i{^)\ickhi)Ag{cihi), 

k=l 1=1 

where the values of A^ that have been already evaluated for the quadrature of 
the univariate integral are reused. For triple and further integrals the same pro- 
cedure can be applied: replace Ag by the Lagrange polynomial of sufficiently 
high degree and then use the resulting analytic expressions for the integrals as 
approximating formulae. Note that also the value of q in (19) is computed by 
Gauss-Legendre with z/ nodes, and thus the same function evaluations of q{x) 
are needed as to compute the different Aq{cihi) in the Lagrange polynomial. 

An alternative way to apply the Filon-type rule is by approximating q{x) 
(piecewisely) by a series over shifted Legendre polynomials (as in done in the 
description of PPM [28,32]): 

u-l 

q{x)^Y.QshtP:iS/h,), 5 = x-x,.i (31) 

By the method of least squares the expressions for the coefficients Qs are 
obtained: 



( 2g + l) rh 



Qs = ^ , ,+1 /„ ?(a^^-i + 5)P:{5/K)d5, m = 0,1,2,.... (32) 



It can then be noted that q = Qo and Ag{6) ^ £a,(5) = - Es=i QshlP*{5 /hi). 
Writing 2>/^^{5) in this form is fully equivalent as using (30), but allows to 
obtain shorter expressions for the formulae approximating the integrals. To 
compute the integrals (32) Gauss-Legendre is used, requiring u function eval- 
uations of q. Suppose we truncate all but the first integral, resulting in a 
method of order four. We need to discretise the integral consistently with the 
order of the method. To this end, we take v = 2. With 

l=i{Z2h), m = VoiZ2h), Z2h = ^Zh = ^{q-X)hi 
and Qs = hl^^Qs, s = l,...,z/ — 1, we then obtain the following 



' f''AMSVoiZ.s)dSr^^^^-'-^^'^^^ 



hjo '' ' '"' ^"' AZh 

'^ Ag{5) (1 + az^s)) d5 ^ r^ Agi5)aZ2s)d5 
Jo 



/ ^Ag{S){l-^{Z2s))dS^- f ^Ag{S)C{Z25)dS^Qif]o + 
Jo Jo 



2Zu 
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which allows us to approximate /g ' B{S)dS. Note that the quadrature approx- 
imation of the non-oscillating integral Jq' Aq{6)d6 vanishes: Jq' 2.AqiS)d6 = 
/q '(g — £,q{6))d6 = 0, since both /q ' £,q{S)d6 and qhi are the Gauss-Legendre 
quadrature approximations of /q ' q{xi_i + (5)(i5. 

To construct a method of order eight, we need to take u = 4 and we have to 
include the double integral. To compute the univariate integral we have now 



1 /•'*'w.^. rry ^,, {Qi + SQ2 + 6Q3)Vo {Qs + Qi){^ + I) + Q2{^ - I) 



- f ^ Aq{6)6r]o{Z25)d6 
hi Jo 



IZh 4Zh 

34(1-0-154(^ + 1) 15Q3^o 
4Z^ ^ + 2Zl 

"' \imZ2s)d6^{Q, + 4 + Qs)m + (^^^ + ^^^^)^° 

sQ2ii+i) + iQi + 6Q3m-i) , 154(1-0 



2Zh 2Zf^ 



Suppose we construct a Magnus method, we consider then the approximation 
of o"2. A similar procedure can be followed to compute the double integral in 
a Neumann method. As in [34] we write the double integral in o"2 as 



'^^ r[B{62),B{6i)]d62d6i = 2 f'^^ ['' Aq{6,)Aq{S2)K,{6i,62)dS2dS,Ui 
Jo Jo Jo 

hi rSi 



+ 2 / ' f " Aq{6,)Aq{62)K2{6i,62)d62d6,U2 

Jo Jo 

+ 2 r ['" Aq{6i)Aq{62)K,{6,,62)dS2dS,U- 

Jo Jo 



4(A-g) 


u 





1 



where Ki{x,y) = yr]o{Z2y) - xr]o{Z2x) , K2{x,y) = ^{Z2x) - i{Z2y), K2.{x,y) 
{x - y)r]o{Z2(x^y)) and 



u, = I ^^^-^^ I , f/. = I ° 2("-^") 

4(A-g-)/ ^ ^ 



The three integrals in (33) are replaced by quadrature by again approximating 
Aq by the polynomial £a • The expression for the third integral is then for 
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instance: 

Jo 

QI-QI-QI-2Q3Q1 -Ql + 15QI - 66QI - A2Q3Qi 

ggi - 405gi - soggQi 225gi\ , /g? - sg^ + egi + 7^4 



+ 4Z3 ^zir^K AZl 

sogggi - 9gi + losgj 225g|\^_ 42g^ + 7ogf + sogj _ s^gi 

4^3 ^ 4Z^ y^ 840Z;, AZl ' 

In practice, one should use a truncated series expansion for small Z^ values 
(see [34]). 

Even higher order algorithms can be constructed including more Magnus (or 
Neumann) terms in the scheme. In [35] a Magnus scheme of order 10 is de- 
scribed where z/ = 5 and in [32] PPM-schemes up to order 16 are presented. 
Note that only the terms where the degree in h is smaller or equal to the re- 
quired degree of the method have to be included in the algorithm, for instance 
in the approximation of a^ and o"4 the term in Q\ can be disregarded in the 
Magnus scheme of order 10. 

The modified Magnus methods and modified Neumann methods are well suited 
for the repeated solution of the initial value problems which appear in the 
shooting procedure. These initial value problems are solved for a fixed po- 
tential q but for different values of A. As shown in [28,34], an A-independent 
mesh can be computed which is then (re) used in all eigenvalue computations. 
Moreover also the value q and the coefficients Qs are computed and stored 
once for all before the start of the shooting process. Algorithm 1 shows the 
basic shooting procedure in which a modified Magnus algorithm is used to 
propagate the left-hand and right-hand solutions. 



5 Some notes on singular problems 



When the problem is singular, either because (a, h) is an infinite interval or 
because at least one of the coefficients p~^, g, w is not integrable up to one 
of the endpoints, then an interval truncation procedure must be adopted. 
Different algorithms are implemented in the available SLP library codes to 
determine a truncated endpoint and appropriate boundary conditions to give 
a prescribed accuracy (see [5]). The SLEDGE package [15] even has algorithms 
for automatically classifying the nature of the problem, regular or singular, 
limit-circle singularity or limit-point singularity and so on. This classification 
information is important to determine whether or not there is a continuous 
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Algorithm 1 A Sturm-Liouville solver based on a modified Magnus method 
1: Use stepsize selection algorithm to construct mesh a = Xo < Xi < ... < 

Xn = h 
2: for i = 1 to n do 
3: Compute q and Qs, s = 1, . . . , z/ — 1 for the ith. interval (Gauss-Legendre 

with V nodes). 
4: end for 

5: Choose a meshpoint x^. (0 < m < n) as the matching point. 
6: Set up initial values for y^ satisfying the BC at a and initial values for y/j 
satisfying the BC at h. Choose a trial value for A. 
repeat 

for ? = to ?T2 — 1 do 

ydxi+i) = e^^'^e'^^^^'^yLixi) 
end for 
for i = n down to m + 1 do 

YRixi^i) = e-'^(^>)e^^>^yij(xi) 
end for 

Adjust A by comparing yiixm) with ynixm). 
until A sufficiently accurate 

spectrum, when there are eigenvalues and how many, and what boundary 
condition should be imposed at a singular endpoint. 

As mentioned before, the algorithm applied in the SLEDGE package relies on 
Pruess coefficient approximation by piecewise constants (namely the midpoint 
values of each interval), and uses repeated extrapolation to achieve accuracy. 
In a ffist pass a crude initial mesh is chosen by an equidistribution process. 
SLEDGE then repeatedly bisects this initial mesh and uses iterated extrap- 
olation. An infinite endpoint is transformed to zero by the (local) change 
of variable t = 1/x, and subsequent bisections near the endpoint are done in 
terms of the variable t. SLEDGE's approach automatically regularizes singular 
endpoints: evaluating the coefficients at the mesh midpoints can be regarded 
as truncating the interval at the midpoints of the ffist and final intervals of the 
mesh. Every time the mesh is bisected these implicit truncation points move 
closer to the singular endpoints. The boundary conditions are always applied 
in the original endpoints. 

For the higher order coefficient approximation methods, a similar approach as 
in SLEDGE can be used to truncate a singular problem. For these methods 
the coefficients are only evaluated in the Legendre nodes. Since the first and 
last Legendre node in an interval are not equal to the beginpoint or endpoint 
of that interval, a singular problem is implicitly truncated. By decreasing the 
size of the first interval (if a is a singular endpoint) or the last interval (if b is 
a singular endpoint) in the mesh, the implicit truncation points move closer 
to the singular endpoint. However no software based on higher order modified 
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integral methods is currently availabe which can automatically handle singular 
endpoints. For this reason, SLEDGE is still the package of choice for many 
applications. 



6 Some experiments 



We consider two well-known test problems. The Coffey-Evans equation is one 
of the more difficult test problems in the literature (test problem 7 in [5], 
introduced in [36]). It is a Schrodinger equation with q{x) = —2j3cos{2x) + 
/9^sin^(2a;) and y{—7i/2) = y{7i/2) = as boundary conditions. The first 50 
eigenvalues for /5 = 30 have been determined. For this potential Aq is close to 
zero and there are very close eigenvalue triplets {A2, A3, A4}, {Ag, A7, Ag}, . . . 
as [3 increases. The second problem is a problem from chemical physics: the 
Woods-Saxon problem [37] defined by 



1 



5t 



3(l+t) 



qix) = -50^^ (34) 

with t = e*^^^^-*/*^'^ over the interval [0,15]. The eigenvalue spectrum of this 
Woods-Saxon problem contains 14 eigenenergies Aq, ..., A13. 

Tables 1-4 show some results for the two test problems obtained with five dif- 
ferent coefficient approximation methods. The first method used is the Pruess 
method of order P = 2. This method is compared with some higher order 
methods which were discussed in section 4.3: a Neumann method and a Mag- 
nus method of order P = 4, and a Neumann method and Magnus method 
of order P = 8. We present for each problem, a selection of the considered 
exact eigenvalues Afc, and the (absolute) error for the corresponding eigenval- 
ues calculated with the coefficient approximation methods. For the moment 
equidistant meshes are used in order to allow easier comparison between the 
different algorithms. An automatic stepsize selection algorithm will be dis- 
cussed afterwards. 

The five different methods were first applied on the same mesh, with 128 steps 
for the Coffey-Evans problem and 64 steps for the Woods-Saxon problem. 
Results are shown in tables 1 and 2. All methods allow the approximation 
of higher eigenvalues or large batches of eigenvalues. However it is clear that 
the higher order methods need much less mesh intervals to reach a prescribed 
accuracy. This can also been seen from the tables 3 and 4, where the number 
of intervals in the equidistant mesh (nint) and function evaluations (nfev) 
are shown that each method needs to reach an accuracy of (approximately) 
10^*^. The data reported in the four tables enable several conclusions: 



17 



Table 1 

Absolute value of (absolute) errors AA^ in the eigenvalues computed for the Coffey- 
Evans problem with different coefficient approximation methods and n = 128 steps 
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le method. 
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Each of the five methods can reproduce accurate results, even for high eigen- 
values. 

However when Neumann or Magnus terms are introduced, less intervals are 
needed in the mesh to reach a certain input accuracy. 
As a result the number of function evaluations is also decreasing with in- 
creasing order of the method. The second order Pruess method only needs 
one function evaluation per interval, while the fourth order methods need 
two and the eighth order methods need four. However to reach the same 
accuracy, the second order method needs many more intervals which causes 
a very high total number of function evaluations. 

Since only a small number of function evaluations is needed, the construction 
of the mesh takes less time for a higher order method. But also the shooting 
process in which the equation is repeatedly integrated at various values of A 
is considerably faster as a result of the smaller number of intervals needed 
in the mesh. The timings shown in tables 3 and 4 were obtained using 
Matlab-implementations of the different methods. 

The modified Neumann and Magnus methods are particularly well suited 
to compute (large) batches of eigenvalues. A remarkably small number of 



Table 2 

Absolute value of (absolute) errors AA^ in the eigenvalues computed for the Woods- 
Saxon problem with different coefficient approximation methods and n = 64 steps 
in the equidistant mesh. aE-6 means a.fO^^. P is the order of the method. 
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mesh intervals is sufficient to be able to compute accurate approximations 
for the higher eigenvalues. There is also a big speed advantage in the fact 
that the repeatedly asked task of integrating the equation at various trial 
values for an eigenvalue is completely separated from the time-consuming 
process of constructing a mesh where all function evaluations are performed 
and saved for later use. 

The fourth order Magnus and Neumann versions reach the same accuracy. 
The eighth order Neumann method seems to be somewhat more precise 
than its Magnus counterpart. This is a consequence of the finite radius 
of convergence of the Magnus expansion (see [23]) which means that the 
steps in the q ^ \ region can not always be taken as large as for the 
Neumann expansion when Zh is large and positive. Note that cases with 
very large positive Z rarely appear in practice. They are e.g. ruled out by 
WKB arguments which usually shrink the interval [38]. 
The Neumann methods are somewhat faster than their Magnus counterpart. 
The evaluation of the matrix exponential in the Magnus method requires 
little extra time for the second order Sturm-Liouville problem. For problems 
with higher dimension, the computation of the matrix exponential may be 
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Table 3 

Absolute value of (absolute) errors AA^ in the eigenvalues computed for the 
Coffey-Evans problem. Different coefficient approximation methods were used on 
an equidistant mesh with nint steps to reach an accuracy of approximately 10~®. 
nfev is the number of evaluations of the potential function and T{s) is the CPU- 
time needed to compute the first 51 eigenvalues. 
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fairly expensive. 
• However truncated Neumann expansions do not respect Lie-group structure 
which can be a limitation in some applications. But as remarked in [21], 
the basic step underlying the Neumann expansion discussed here is the 
transformation (20), which always preserves Lie-group structure. Departures 
from a Lie group might occur only in the function u, in other words in the 
correction term. This results in far less severe loss of Lie-group structure 
than is the case with classical Runge-Kutta or multistep methods. 

Of course using a uniform mesh is rarely a good idea e.g. when dealing with 
(truncated) singular problems. For automatic software a stepsize selection al- 
gorithm should be used. We refer to [15] for the procedure used in SLEDGE. 
A meshing algorithm has also been proposed for the PPM in e.g. [28] and for 
the numerical solution of regular Schrodinger problems with an eighth order 
Magnus method in [34]. We briefly describe a procedure which can be used 
for the eighth order Neumann method and apply it on the Woods-Saxon prob- 
lem. The resulting mesh is shown in Figure 1, and the results obtained over 
this mesh are listed in Table 5. 

The stepsize selection algorithm applied here is based on a local error estimate. 
Let tol be an input tolerance parameter. In order to attain the local error es- 
timate ej = tol over the interval i, the stepsize hi is chosen as a function of 
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Table 4 

Absolute value of (absolute) errors AA^ in the eigenvalues computed for the 
Woods-Saxon problem. Different coefficient approximation methods were used on an 
equidistant mesh with nint steps to reach an accuracy of approximately 10^^. nfev 
is the number of evaluations of the potential function and T(s) is the CPU-time 
needed to compute the 14 eigenvalues. 
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Fig. 1. The mesh resulting from the stepsize selection algorithm of the eighth order 
Neumann method for the Woods-Saxon potential with input tolerance tol = W~^: 
the values of the potential q{x) at the mesh points are marked by dots. 

the previous stepsize /ij_i as follows. First we compute 



hi = hi — j , h = hi^i 



(35) 
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Table 5 

Absolute value of (absolute) errors AA^ in the eigenvalues computed for the Woods- 
Saxon problem. An eighth order Neumann method is used on the mesh (with 46 
steps) shown in figure 1. 
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where e is the error estimate. A decision is taken in terms of A = \hi/h — 1|. 
If A > 0.1 the procedure is repeated with /i = /ij. If A < 0.1, h is accepted as 
the stepsize for the interval i. 

To construct the error estimate e we consider the maximal (absolute) difference 
between the entries of the transfer matrix of an (embedded) sixth order scheme 
and the entries of the transfer matrix of our eighth order algorithm. This 
difference is evaluated through the sum of all terms in Q^. When scanning 
over A we have used only three values, those such that Z(hi) = {q — A)/if = 
— m^vr^, m = 0,1, 2. The selection of only these was mainly intended to speed 
up the evaluation but this is enough as the error decreases like 0{l/\/X) for 
large Z- values [26,28], and as confirmed by experimental tests showing that 
the error is indeed larger for smaller values of Z. 

Let us consider now a singular problem, just to illustrate that the modified 
integral series methods can be extended to singular problems. The general 
Woods-Saxon problem is a Schrodinger problem of the form 



y"{x) 



'l{l + l] 



x^ 



+ lix) - A y{x) 



with q{x) as in (34) and x G [0, -|-oo]. When the orbital quantum number 
/ equals zero, the potential is a well behaved, nonsingular function, as the 
Woods-Saxon problem we considered before. When / > the problem is sin- 
gular in the origin. We will compute the eigenvalues for the problem with 
/ = 2. This problem has 13 eigenvalues in its spectrum. The infinite endpoint 
can be dealt with by a change of variable converting (implicitly) to a finite 
interval (as in SLEDGE) or with explicit interval truncation using WKB ar- 
guments as in [38]. We describe here only a way to deal with the singularity 
in the origin and consider the (truncated) problem over the interval [0,20]. 
We apphed the stepsize selection algorithm discussed above over the interval 
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Table 6 

Absolute value of (absolute) errors AA^ in the eigenvalues computed for the singular 
Woods-Saxon problem. An eighth order Neumann method is used with user input 
tolerance tol = 10"'^. nint is the number of intervals in the initial mesh and T is 
the CPU-time needed to compute all eigenvalues. 
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[e, 20] and added to the resulting mesh the interval [0, e]. Note that the poten- 
tial function is only evaluated in the Legendre points of the interval [0, e] and 
not in the singular endpoint a = 0. A first approximation of an eigenvalue is 
computed over the mesh. Then the interval [0, e] is bisected and a further eigen- 
value approximation is computed. The process of bisection in [0, e] is repeated 
until two successive eigenvalue approximations agree within the user specified 
tolerance. At each iteration, the shooting algorithm for the next eigenvalue ap- 
proximation is started using the approximation last obtained. Table 6 shows 
the results for two different values of e. A larger value of e requires of course 
more bisections (nbisec) of the interval [0,e]. The results obtained are within 
the requested accuracy tol = 10"'''. 



7 Conclusion 



In this paper we discussed some techniques which allow the efficient approxi- 
mation of high eigenvalues of a Sturm-Liouville problem. We focused in partic- 
ular on algorithms based on approximation of the coefficient functions of the 
differential equation. The simplest coefficient approximation method is the 
so-called Pruess method, which replaces the coefficient functions over each 
mesh interval by the midpoint value in that interval and then solves (analyt- 
ically) the approximating problem. This Pruess method has some significant 
advantages over shooting methods based on standard initial-value solvers es- 
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pecially when looking for higher eigenvalues and was implemented in the well- 
known Sturm-Liouville solver SLEDGE. The oscillations in the eigenfunctions 
no longer determine (restrict) the step sizes. Now the step sizes depend on 
the errors made in replacing the coefficient functions by piecewise constant 
approximations. It is clear that larger steps could be taken when the coeffi- 
cient functions are replaced by higher order polynomials. However when the 
coefficient functions are replaced by polynomials of degree greater than zero, 
the approximating problem is not really easier to solve than the original prob- 
lem. Therefore only piecewise constant (and linear) polynomial approxima- 
tions were used for a long time. But, Neumann or Magnus integral series offer 
a way to construct methods based on higher order (piecewise) polynomial ap- 
proximation. Using (modified) Neumann or Magnus schemes, we can construct 
methods which still allow the easy analytic integration of an approximating 
problem with piecewise constant coefficients but use higher order polynomial 
approximations to construct some extra (correction) terms. Depending on the 
number of terms included in the Neumann or Magnus series, algorithms of 
different orders can be constructed. Experiments show that indeed these Neu- 
mann and Magnus integrators share the advantages of the Pruess method and 
allow to approximate high eigenvalues in a remarkably small number of steps. 
Also a singular problem was considered to illustrate that singular endpoints 
can be dealt with in a relatively simple way. 
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